pro cels_plot_optimized, sim_name

  StartY=cels_get_parameter(sim_name, 'STARTY')
  EndY=cels_get_parameter(sim_name, 'ENDY')
  pollutant=cels_get_parameter(sim_name, 'POLLUTANT')
  stations=cels_get_parameter(sim_name, 'STATIONS')
  stations_all=[stations, cels_get_parameter(sim_name, 'BG_STATIONS')]
  restore, cels_filestr(/Input, sim_name + '/state.sav')
  restore, cels_filestr(/Input, sim_name + '/regions.sav')
  restore, cels_filestr(/Input, sim_name + '/sensitivity.sav')
  restore, cels_filestr(/Input, sim_name + '/optimized.sav')

  nStations_all=n_elements(stations_all)

  color=['BLU6', 'YGB4', 'PUR5', 'PUR7', 'RED7']
  colorBG=['BLU2', 'YGB2', 'PUR2', 'PUR2', 'RED2']

  date_label=LABEL_DATE(Date_format=['%Y'])

  !p.multi=[0, 2, 3]
  position=plot_position(2, 3)

  ;Get measurement time inverval
  y_temp=y_time[where(y_station eq 0)]
  dt=median(y_temp[1:-1]-y_temp[0:-2])
  undefine, y_temp

  ;plot optimized mole factions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ps_start, filename=cels_filestr(/Output, sim_name + '/mf_optimized.eps'), /quiet
    for si=0, nStations_all-1 do begin
      wh=where(y_station eq si and y_time ge julday(1, 1, StartY, 0))

      yrange=[min([y[wh], y_opt[wh]]),max([y[wh], y_opt[wh]])]
      
      cgPlot, y_time[wh], y_opt[wh], xtickformat='LABEL_DATE', xtickunits='Year', $
        yrange=yrange, thick=3, /nodata, position=position[si, *], charsize=1.8
      
      y_polyfill=[y[wh] - y_error[wh], reverse(y[wh] + y_error[wh])]
      wh_polyfill=where(y_polyfill ge yrange[1])
      if wh_polyfill[0] ne -1 then y_polyfill[wh_polyfill]=yrange[1] - 0.01*(yrange[1] - yrange[0])
      wh_polyfill=where(y_polyfill le yrange[0])
      if wh_polyfill[0] ne -1 then y_polyfill[wh_polyfill]=yrange[0] + 0.01*(yrange[1] - yrange[0])
      cgColorfill, [y_time[wh], reverse(y_time[wh])], y_polyfill, color=colorBG[0]

      y_gaps=plot_gaps(y_time[wh], y[wh])
      cgPlot, y_gaps[0], y_gaps[1], color=color[0], /overplot, thick=4
      y_gaps=plot_gaps(y_time[wh], y_opt[wh])
      cgPlot, y_gaps[0], y_gaps[1], color=color[4], /overplot, thick=4

      cgText, position[si, 0] + 0.02, position[si, 3] - 0.05, Stations_all[si], /normal, charthick=3
      
    endfor
    
    cgText, 0.02, 0.5, strcompress(pollutant + ' (pmol mol!u-1!n)'), /normal, charthick=3, orientation=90, alignment=0.5
  ps_end



  for lri=0, N_LR-1 do begin

    restore, cels_filestr(/Input, sim_name + '/ar_' + string(lri, format='(I02)') + '.sav')
    restore, cels_filestr(/Input, sim_name + '/emissions/lr_' + string(lri, format='(I02)') + '_2008.sav') ;MLR
    
    lon_domain_lr=reform(lon_domain[*, lri])
    lat_domain_lr=reform(lat_domain[*, lri])

    area=areagrid(lonq, latq)

    q=q*365.25*24.*3600.*area  ;kg/year
    
;    q_opt=q
    q_opt=fltarr(EndY - StartY, n_elements(lonq), n_elements(latq))
    q_ap=fltarr(EndY - StartY, n_elements(lonq), n_elements(latq))
    q_opt_mean=fltarr(n_elements(lonq), n_elements(latq))
    q_ap_mean=fltarr(n_elements(lonq), n_elements(latq))

    caldat, x_time, dummy, dummy, x_year

    for yi=StartY, EndY-1 do begin
      x_temp=x[where(x_lr eq lri and x_year eq yi)]
      x_ap_temp=x_ap[where(x_lr eq lri and x_year eq yi)]
      q_ap[yi - StartY, *, *]=q ;CHANGE FOR IAV AP!!!!!!!!!!!!!!!!!!! MLR

      q_opt_temp=fltarr(n_elements(lonq), n_elements(latq))
      for xi=0, n_elements(x_temp)-1 do begin
        q_opt_temp[ar_wh[xi]]=q[ar_wh[xi]]*x_temp[xi]/x_ap_temp[xi]
      endfor
      q_opt[yi - StartY, *, *]=q_opt_temp
      q_opt_mean+=q_opt_temp
      q_ap_mean+=q  ;CHANGE FOR IAV AP!!!!!!!!!!!!!!!!!!! MLR
    endfor

    q_opt_mean=q_opt_mean/(EndY-StartY)
    q_ap_mean=q_ap_mean/(EndY-StartY)
    
    q_diff=q_opt_mean-q_ap_mean
    
    lonrange=[min(lon_domain_lr), max(lon_domain_lr)]
    latrange=[min(lat_domain_lr), max(lat_domain_lr)]
    if lonrange[1] ge 180. then lonrange[1]=179.9
    
    diff_levels=fltarr(11)
    diff_max=max(abs(q_diff), /nan)
    diff_max_order=fix(alog10(diff_max))
    
    diff_levels[6:10]=10.^(findgen(5) - 5. + diff_max_order + 1.)
    diff_levels[0:4]=reverse(-1.*10.^(findgen(5) - 5. + diff_max_order + 1.))
    diff_labels=[string(reverse(findgen(5) - 5. + diff_max_order + 1.), format='(I2)'), '+!c-!c', $
      string(findgen(5) - 5. + diff_max_order + 1., format='(I2)')]

;    q_diff[where(q_diff eq 0.)]=!values.f_nan

    mr_plot_map, q_diff, lonq, latq, filename=cels_filestr(/Output, sim_name + '/' + string(lri, format='(I02)') + '_q_diff.eps'), $
      /square, limit=[latrange[0],lonrange[0],latrange[1],lonrange[1]], /diff, ct=22, /reverse_ct, $
      cbtitle='Optimized - prior  (Log!d10!n(kg / year))', levels=diff_levels, labels=diff_labels

    mr_plot_map, alog10(q_ap_mean), lonq, latq, filename=cels_filestr(/Output, sim_name + '/' + string(lri, format='(I02)') + '_q_ap.eps'), $
      limit=[latrange[0],lonrange[0],latrange[1],lonrange[1]], ct=25, /reverse_ct, $
      cbtitle='Prior emislrions  (Log!d10!n(kg / year))' ;, range=range

    mr_plot_map, alog10(q_opt_mean), lonq, latq, filename=cels_filestr(/Output, sim_name + '/' + string(lri, format='(I02)') + '_q_opt.eps'), $
      /square, limit=[latrange[0],lonrange[0],latrange[1],lonrange[1]], ct=25, /reverse_ct, $
      cbtitle='Optimized emissions  (Log!d10!n(kg / year))' ;, range=range
    
  endfor


end